function is_ice = mycode_is_ice(flag_rcp, flag_model, flag_period)

%%% Distinguish between ice-free and ice-covered lakes
%%% We define the ice-covered lakes as those with more than half of the ensemble members showing that the multi-year averaged ice duration 
%%% is longer than 1 day·year-1 or longer than 1 day·season-1 at the annual or seasonal scale

% flag_rcp: The difference between the flag_rcp (pre-industrial or future projections) and present-day period 
% #2 = picontrol        #3 = rcp26      #4 = rcp60      #5 = rcp85

% flag_model: The results of specific lake model or ensemble mean
% #1 = albm     #2 = lake      #3 = simstrat-uog		#4 = vic-lake		#5 = ensemble mean

% flag_period: The seasonal or annual results
% #1 = DJF       #2 = MAM       #3 = JJA      #4 = SON      #5 = annual


ice_duration_rcpMean_ensemble = evalin('base','ice_duration_rcpMean_ensemble');
rcp = evalin('base','rcp');
period = evalin('base','period');

flag_var = 1;   % 'ice_duration_mean'

if flag_model < 5
    
    present_ice_duration = ice_duration_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ flag_model }{ flag_var };
    future_ice_duration = ice_duration_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ flag_model }{ flag_var };

else
    
    present_ice_duration = cat(3, ice_duration_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 1 }{ flag_var }, ...
                    ice_duration_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 2 }{ flag_var }, ...
                    ice_duration_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 3 }{ flag_var } );
    future_ice_duration = cat(3, ice_duration_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 1 }{ flag_var }, ...
                    ice_duration_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 2 }{ flag_var }, ...
                    ice_duration_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 3 }{ flag_var } );
end

if flag_rcp == 2
    a = future_ice_duration;
    future_ice_duration = present_ice_duration;
    present_ice_duration = a;
else
end
clear a

is_ice = (present_ice_duration + future_ice_duration)/2;
is_ice(is_ice<=1) = 0;
is_ice( is_ice~=0 & ~isnan(is_ice) ) = 1;
is_ice = mode( is_ice,  3 );
is_ice(is_ice==0) = nan;

end




